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Vh ' Abstract 

The photon density operator function is used to calculate light 
beam propagation through turbulent atmosphere. A kinetic equation 
for the photon distribution function is derived and solved using the 
method of characteristics. Optical wave correlations are described in 
terms of photon trajectories that depend on fluctuations of the refrac- 
tive index. It is shown that both linear and quadratic disturbances 
produce sizable effects for long-distance propagation. The quadratic 
terms are shown to suppress the correlation of waves with different 
wave vectors. 
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We examine the intensity fluctuations of partially coherent beams 
(beams whose initial spatial coherence is partially destroyed). Our cal- 
culations show that it is possible to significantly reduce the intensity 
fluctuations by using a partially coherent beam. The physical mech- 
anism responsible for this pronounced reduction is similar to that of 
the Hanbury-Braun, Twiss effect. 



1 Introduction 

The study of the interaction of light beams with a random media is of great 
importance for applications in such areas as astronomy, laser communica- 
tion, laser radar systems, etc. Fluctuations of the atmospheric refractive 
index, caused by turbulent eddies, affect electromagnetic waves. The effects 
of turbulence have been investigated both theoretically and experimentally 
over the past few decades. (See, for example, the monographs [I]-j3] and the 
surveys |l]-[7j.) Additional beam spreading, beam wandering ("dancing"), 
intensity fluctuations (known as "scintillations"), etc. caused by turbulence 
limit significantly the range and the performance of free-space communica- 
tion systems. 

It has been established that coherent laser beams are very sensitive to 
atmospheric conditions. Realizing this, the idea of utilizing of partially co- 
herent beams (PCB) for practical purposes has arisen. It was shown by 
numerous researchers (see, for example, publications |H]-C1]) that partially 
coherent beams are less affected by turbulence than fully coherent beams. A 
specific case (in which the spatial coherence of the signal-carrying laser beam 
is partially destroyed before it is launched into the atmospheric channel) was 
considered by many researchers. This beam has the angular spread in free 
space of the order of \/lc (A and Ic are the wavelength and the transverse 
correlation length of phase distortion at the source aperture, respectively), 
which is larger than that of the coherent beam. For not too large distances 
of propagation, beam spreading due to the small initial coherence length, 
/c, may dominate throughout the trajectory. One can say that the effects of 
turbulence are masked by the larger initial free-space spreading. At the same 
time it is evident that, with increasing propagation distance, the atmospheric 
inhomogeneity becomes the dominating factor. This results in a beam size 
that is almost independent of the initial correlation length. 

The dependence of the intensity fluctuations on the initial coherence is 



more intricate than the dependence of beam spreading. Provided the de- 
pendence of the intensity fluctuations on the initial coherence is given, the 
noise/signal ratio can be controlled by choosing the optimal initial coherence 
length, Ic- This tempting opportunity has stimulated many studies in this 
field P, [ini-EIl- It was shown that control is possible for small fluctuations 
(weak turbulence or short-distance propagation). 

At moderate and strong fluctuations the situation is more complicated. 
Banakh et al. jTH] have shown that for long distances of propagation or for 
strong turbulence the intensity fluctuations have a tendency to saturate at a 
level which depends on the initial spatial coherence. This level may be much 
lower than that of the coherent beam by the factor, (ro/lc)'^, where tq is the 
initial radius of the beam. The paper deals with the physical model when 
the measuring device (detector) has the response time, r^, greater than the 
coherence time, Tg, of the radiation at the beginning of the trajectory. In 
this case, the detector averages the signal over the initial fluctuations, which 
are due to temporal intensity fluctuations of the source and (or) may be gen- 
erated when a coherent beam is transformed into PCB. The authors of ^Hl 
have solved an equation for the fourth-order coherence function (the fourth 
moment of the field). The derivation of the equation for fourth moment in 
the limit of Markov approximation is given in P], ^. Its solution is based 
on the approximate method developed by Yakushkin in [221 • The authors of 
[TH] have modified the Yakushkin approach to be applicable to the case of 
PCB. The Yakushkin method is based on the observation that, for long dis- 
tance propagation, the dominant contribution to the fourth-order correlation 
function comes from the products of two second-order coherence functions 
just as if Gaussian statistics for radiation field were valid. The Dashen anal- 
ysis [23j in terms of the path integral formalism and the results of Fante j21] 
obtained by employing the extended Huygens-Fresnel principle have justified 
the Yakushkin idea. It seems to be quite reasonable to consider that after 
long-distance propagation through a medium with a random refractive index, 
the statistics of the photon flux acquires some of the properties of thermal 
radiation. 

For similar conditions, when the coherence time of the quasimono chro- 
matic laser source is much smaller than the detector's integration time inter- 
val, Td, Korotkova et al. J2I] have derived an analytical expression (see Eq. 
(33) in J2I]) for the intensity fluctuations. The paper [21] has generalized the 
analytical approach developed by Andrews and Phillips (see monographs [2] 
and pi) to the case of PCB. This approach is a modification of the Rytov 



theory. Namely, in addition to the first-order perturbation of the complex 
phase of the field caused by refractive index fluctuations (as in Rytov the- 
ory), second-order perturbation terms are taken into account. The results of 
PT] differs from those in ^B]- More comments are presented in Section 4. 

The purpose of our paper is to develop a new quantum perturbation 
theory based on path integrals. This approach allows us to obtain the scin- 
tillation index for the case of a PCB, including a limit of a strong turbulence 
and long-distance propagation. Our approach uses the distribution function 
of photons (for photon density in phase space), /(r, q, t), where r and q are 
the coordinate and the momentum of photon, to describe the beam char- 
acteristics at an arbitrary instant, t. The first and second moments of the 
distribution function (< / > and < // >) are used to obtain the beam size 
and intensity fluctuations, respectively. Also, higher-order moments will be 
obtained for long propagation paths. 

The next Section deals with the definition of the distribution function. 
The equation governing the evolution of the distribution function will be 
derived. 

2 The photon distribution function and its 
evolution 

The Hamiltonian of photons in a medium with a fluctuating refractive index 
is given by 

H = J2 ^^k&k &k - Yl ^^k^k'^k &k+k' , (1) 

k k,k' 

where the two terms on the right-hand side describe photons in a vacuum and 
the effect of refractive index fluctuations, respectively; b^ and 6k are creation 
and annihilation operators of photons with momentum k, ^cuk = hck is the 
photon energy, c is the speed of light in a vacuum, and rik is the Fourier 
transform of the refractive index fluctuations 6n{r) . The Fourier transform 
is defined by 

n^ = ^jdVe''''5n{v), (2) 

were V = L^LyLz is the normalizing volume. 

Eq. (^ follows from the representation of the energy of the random 
medium + electromagnetic field given in J2S] (Chapter 15) in the limit of 
small wave-vectors k' {k' <^ k) and of atmosphere refractive index close to 



unity (^^(r) — 1 << 1). The former means that the scale of spatial inhomo- 
geneity of turbulence is much greater than the wavelength of the radiation. 
For simplicity, we consider here only the case of polarized light with a fixed 
polarization throughout the distance of propagation. Depolarization effects 
due to atmosphere turbulence are very small. (See, for example, papers [211 
and 27 .) Also, the terms describing the zero-point electromagnetic energy 
are omitted in Eq. (P). 

The photon distribution function is defined by analogy with the distribu- 
tion functions for electrons, phonons |28j, etc. It is given by 

/(r,q,t) = ^Ee"^'^&WVk/2, (3) 

k 

The distribution function will be used to describe beams with character- 
istic sizes much larger than the photon wave-length. Thus, we will restrict 
a sum over k by some ko {k < ko << Qq, where go is the wave vector corre- 
sponding to the central frequency Uq of radiation, uq = cqo). At the same 
time ko is chosen to be large enough to sample the spatial variation of the 
light intensity. 

After integration of the operator function /(r,q, t) over the volume V, 
we obtain the operator for the total number of photons with momentum q: 

ldVf{r,ci,t) = b^b^. (4) 

Similarly, the quantity obtained after a summation of /(r, q, t) over q may 
serve as a photon density averaged over a small spatial area with the size vr/fco- 
This is very similar to the Mandel operator [22] (Chapter 12) introduced in 

Here and in the remainder of this paper, we use the Heisenberg represen- 
tation for all operators. Thus, the evolution of /(r, q, t) is determined by the 
commutator with the total Hamiltonian: 

5,/r,q,t) = ^[/r,q,t),i7]. (5) 



Eq. (0) can be written explicitly as 
<9t/(r, q, t)+Cq9r/(r, q, t)-iuJo JZ e'^^^^n^ 

k 



/(r,q+|,t)-/(r,q-^,t' 



= 0, 
(6) 



where Cq = SqtUq. 

Considering the characteristic values of the photon momentum to be much 
greater than the wave vectors of turbulence, we can express the difference 
of functions in square brackets by the corresponding derivative. A detailed 
discussion of this approximation is given in Section 6. Then, after summing 
over k, we obtain 

{9i + Cq9r + F(r)aq}/(r,q,t) = 0, (7) 

where F(r) = u!odrn{r). As we see, the photon distribution function is gov- 
erned by a kinetic equation in which the random force F(r) originates from 
atmospheric turbulence. 

The distribution function determines density of photons at a point (r,q) 
of the phase space at time t. Eq. (jTj) may be interpreted as the equation 
governing the evolution of a particle distribution function in which the state 
of each particle is described by the individual coordinate r and the momen- 
tum q. The trajectories of these particles may be obtained from the solution 
of the equations of motion: 

dr(t) , , ,, 



dt 

dqit) 



F(r(t)). (8) 



dt 
Then, the general solution of Eq. ((7j) is given by 

where the function 0(r, q) is the "initial" value of /(r, q, t), i.e. 

''^ k k 

and the "trajectories" r(t') and q(t') pass through the point r, q at t' = t 
[i.e. r(t' = t) = r, q(t' = t) = q)]. As one can see, the photon distribution 
function at an arbitrary instant t is expressed via the operators b't, 6q defined 
for some fixed to (^o is chosen to be equal to in Eq. Q. It is convenient to 
put t — to = z/c. Thus, to is the instant when photons exit from the source. 



The initial photon statistics (at to); determined by the source properties, is 
assumed to be given. 

We consider here the propagation of hght beams with narrow spread 
(paraxial beams). In this case, k±,q± -C go, where index (_l) means perpen- 
dicular to the direction of propagation (the z-axis) components. The relative 
effect of turbulence on q^ is negligible because of the large value of go- At the 
same time, q±, which determines a beam divergence, can be increased con- 
siderably due to turbulence (compared to the initial value). Therefore, beam 
characteristics should be modified significantly for the case of long distance 
propagation. 

It follows from Eq. (jS)) that the evolution of transverse photon momentum 
is given by 

ci±{t') = q^ + l'dt"FMt")]. (11) 

Similarly, we obtain an expression for r(t') 



r(t') = r - Cq(t -t')-- f dt"{t" - t'WMt")]- (12) 

Qn Jt 



qoJt 
Then, Eq. (jH)) can be written as 



f{v,q,t) = Jr-c^t + - fdt't'F4r{t%q- f dt'F^rit')]]. (13) 

t go "'0 Jo J 

A regular iterative procedure is applicable here to expand r(t') in powers 
of F. Then, substituting the explicit terms r(t') into Eq. (fT!^ . we will obtain 
solution of the problem. In particular, the first and second moments of /, 
which describe beam spreading and intensity fluctuation, can be calculated. 
Applying this perturbation method, we can investigate the effects of the 
initial partial (spatial) coherence on beam spreading and scintillations. 

3 Beam spread and intensity fluctuations 

The intensity of radiation in the 2;-direction at r can be presented in the form 

Iir) = cJ2hco^fiv,q,t). (14) 



This can be rewritten as 

(15) 
Let us restrict ouselves to only linear in F terms in Eq. (J12p . In other words, 
we put r — Cq(t — t') for r(t') in arguments of Fj_ of Eq. p5|) . After changing 
variables q_L — /g dt't'F±[r{t')] — >• q_L and using the relation z = ct, Eq. (jT^ 
is transformed to 

q,k 

The stochastic variables F and ^klq)? which are of different nature, are sep- 
arated in Eq. ()16|). Averaging of each factor in the sum can be performed 
independently because of the absence of correlations between the source fluc- 
tuations and the refractive index fluctuations. Thus, we have 

< /(r) >= chuoYl < e-^'^^^^^-'i^(^/^o)+(^/^«)/o^*'*'^^I^-^('i)*'l> >< 0k(q) > • 

q,k 

(17) 
In Eq. fITTj) and throughout this paper, we shall calculate average values of 
functionals of 6n. This can be carried out when the statistics of 6n is known. 
Usually, 6n is assumed to be a Gaussian random variable with known covari- 
ance < 6n{r)6n{r') >. The covariance is defined by its Fourier transform, 
ip{g), with respect to the difference r — r'. The dependence ^(g) is often 
approximated by the von Karman formula 

u \ nnQQr-2ea;p[-(^/o/2vr)2] 

^(g) = omsc^ [g^ + L,Y/^ • ^^^^ 

The quantities Lq and /q are the outer and inner scales sizes of the turbulent 
eddies, respectively. In atmospheric turbulence, Lq may range from 1 to 100 
meters, and /q is usually on the order of several millimeters. C^ is known 
as the index-of-refraction structure constant. In most physically important 
cases the quantity Lq ^ in the denominator of Eq. (fTS|) can be omitted. In 
this case, the von Karman spectrum is reduced to the Tatarskii spectrum P] . 
Using the explicit form for the turbulence fluctuations, Eq. fll7|) becomes 

< /(r) >= chuJo J2 e-^'^^['-^-'i^("/5°)l-'=x-'^ < 0^(q) >, (ig) 

q,k 

8 



where the effect of turbulence is represented by a quantity T {T = 0.558C^/(^ ) 
When obtaining Eq. (fT^ we have assumed that the distance of propagation 
is much greater than the characteristic length of the turbulence. The con- 
dition z » I/Q is sufficient to satisfy this requirement for any regime of 
propagation. Also, the Tatarskii expression for the turbulence spectrum was 
used. 

Now we consider the average value of < 0k (q) >• It depends entirely on 
the source properties. Let us consider a source field with the following mode 
structure 




E'ij) = i: ^-^ [e''"''^nMbn + e-"^"'KMbt], (20) 



where the normalized function, $(r_L), wave vector, g„, and frequency, Un, 
describe the profile and eigenfrequency of n-th mode, respectively. This field 
should be matched with the field in the atmosphere, 

E^'-ir) = Y: (^^) ' '[e'^'^^&c + e-^^%], (21) 

in the plane of the transmitter, where functions $„(rx) are assumed to be 
known. When we deal with single-mode laser radiation (for example, with 
the n = mode), only one term in the sum in Eq. (j2m) should be retained. 
The other terms can be omitted. It is important to point out that rigorous 
matching conditions would also involve the vacuum fields to maintain correct 
commutation relations for all field operators. Nevertheless, the vacuum fields 
may be neglected in our case for two reasons: (i) the photon detectors are 
assumed to be of the absorbing type, hence they are not sensitive to vacuum 
fields, and (ii) we consider only linear (in E) propagation of the radiation. 
In this case, 

V,,o = h{L,Ly)-^'^ j dv^e-'^^'^^{v^), (22) 

where the index in $ is dropped for brevity. 

Until now, all phase distortions were not considered. In practice, stochas- 
tic phase distortions may be introduced by means of a rotating phase diffuser 
placed in front of the aperture. Mathematically, the effect of the phase dif- 
fuser may be taken into account by introducing the multiplier e"*'^^'"^) (see, 
for example, ^Tj) into the integrand of Eq. (|22j) . where '~p{'r±) = ar^ and 



a is a Gaussian random variable with covariance < 



"XjyJ 



>= \-\ Then, 



considering $ to be a Gaussian-type function ($ = -j^e ''"±/'"o)^ we obtain 



< 0k(q) >= 






y^ro 



< b-^b > e-^'(''o/8)-9'(^i/2) 



(23) 



where rf = rl/(l + 2r'^X^ ^) and < 6+6 >= 



for the coherent state |/5 > of 



the laser radiation. Here, symbols q and k are the perpendicular components 
of the wave vectors. As we see, the effect of partial coherence is represented 
by the value of rf. In the limiting case of Ac 
opposite case of small correlation length Ac - 

the quantity b = rf/rl is a measure of a spatial coherence of the beam. 
Using Eqs. (0 and ((221), it is found that 



->■ oo, rf tends to Tq. In the 
0, rf tends to A^/2. Hence, 



< Iiz,v± = 0) >=/o 



1 + 



Az' 



qlrlrl 



+ 



k^T 



(24) 



where /q is equal to < /(r) > at r_L = and 2; = 0. Eq. (j24j) coincides with 
the corresponding Eq. (39) of jlj when C\ does not depend on distance z 
and r\ = vq. 

The intensity for arbitrary r^ can be obtained from (Eq. 1^^ by multi- 
plying its right-hand side by the factor exp{ — ^[1 + i% a + ^^1^^}. The 
average beam radius R defined as. 



R' 






is given by 



i?^ = ^ 



1 + 



Az' 



q^r^rf 



+ 



Iz^T 



(25) 



(26) 



This coincides with Eq. (4) of the paper ^1], where the effect of partial 
coherence was studied. As one can see, only the second term in the square 
brackets depends on the initial coherence via r^. This term describes the 
diffraction spreading of the beam in free space. It depends on both the initial 
beam radius rg/v^ and the coherence length l^ via ri. The spreading may 
be enhanced considerably, if Ac is decreasing. In this way, the diffraction 
divergence may exceed the divergence due to turbulence (the third term) 
for a broad range of distances. In this case, we note the independence of 



10 



(R/^c) 



10^ 



10' 



10 



10^ 



10 




I 1 1 11 1 1 — I I I I 



IQ-' 



— 1 — I I I 1 1 1 1 



10"' 



i'A) 



10' 

2 



10-^ 



10" 



Figure 1: The dependence of beam radius on the coherence length parameter 
('^i/'^o)^ for % = lO'^m"^, (/o/27r) = lO^^m. There is a wide range of (ri/ro)^ 
for which the beam radius is almost constant. 

the beam radius on the turbulence strength, i.e. on the weather conditions. 
Nevertheless, it follows from Eq. (|24p that the "turbulent" term dominates 
in the limit oi z ^ oo. Fig. 1 shows the dependence E?{r\/rQ) for two 
different distances and turbulence strengths. 

Prior to considering the intensity fluctuations, it is useful to analyse qual- 
itatively peculiarities of the wave correlations in the course of their propa- 
gation through the atmosphere. It follows from Eqs. (j23p and (jl9|) that 
the characteristic values of k contributing to < / > are less than or of the 
order of the smallest of the quantities ^^, 2^, [z^T]^^/"^. This means that 

the two waves, &q+k/2 and &q-k/2, are correlated (< ^q+k/2^q-k/2 > 7^0) if 
k satisfies the above requirement. In the limit of 2; — ;> oo, each wave corre- 
lates with itself only as if the beam originates from a thermal source. The 
characteristic distance for wave randomization can be obtained from the re- 
quirement that the "turbulent" term should be dominant term in Eq. ()26|). 
The corresponding criteria are given by 

These inequalities distinguish the range of strong turbulence, which is of the 
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most interest for our studies. 

The intensity fluctuations are determined by the expression 

^ q,kq',k' 

(27) 
where the symbol { : : } means the normal ordering of the creation and anni- 
hilation operators. (See more detail in |29].) 
Let us introduce the notation: 

^(q, q'; k, k') = &q+k/2^q'+k'/2 V-k'/2 Vk/2. 

In the limit z —* oo we have 

< G > = < feq+k/2Vk/2 >< ^q'+k'/2V-k'/2 > (28) 

+ < ^q+k/2V-k72 >< ^q'+k'/2^q-k/2 >= '^q"'q'^k,0'^k',0+'^q+k/2"'q-^k/2^q,q'^k,-k' 

where n^ =< b^b^ >. The terms which describe two pairs of waves with 
coinciding indices in each pair, have nonzero values. At large but finite z, 
correlation of waves with somewhat different indices ( "nondiagonal" terms) 
may also occur. These are terms with: (i) k, k' < (z^T)~^/^ or (ii) |q — q' + 
(k + k')/2|, |q' - q + (k + k')/2| < (z^T)-^/^. The estimates, (i) and (ii), 
follow from the requirement that the deviations of the wave-vectors from 
their "diagonal" values are less than or of the order of the reduced beam 
radius. The last condition is determined by the "turbulent" term in Eq. (j26|) 
for the case of long distance propagation. The intersection of both ranges of 
wave-vectors confined by inequalities (i) and (ii) may be neglected because 
of small volume in wave-vector space. Then, Eq. ()28|) is reduced to 



<: I\r) :>= ^^^^(Si + S2) < G(q, q'; k, k') >, (29) 



where Si^2 means summation, with the restrictions (i) and (ii), respectively. 
The second sum is reduced to the first one by renaming the indices. Then 
considering the quantity (l/^)&q+k/2^q-k/2 as a spatial Fourier component 
of the distribution function, we have 

<• ^^(^^ ■> = S < e-^^^[^-'^('l)*l+'^'['-'^('l')*l + (^/'?o)/Jdt't'(kF[r(q,t')]+k'F[r(q',t')])} 

(30) 
12 



: 0k{q- rrft'F[r(q,t')]}0k'{q'- T rft'F[r(q', t')]} : 
JO Jo 



>, 



where the summation is with restriction (i). For the sake of brevity, we again 
denote by F the component of the force perpendicular to the z-axis. 

The averaging of the product : (jxp '■ with respect to source variables can 
be performed straightforwardly. (See the derivation of Eq. (I23|l . ) It is given 
by 

<: MQ)MQ') ■■>= {yiS'} ^ ^^^^^^ ^ ^^^^ 

Xe-(Q-Q')'^o/4-(Q+Q')'r-2/4-{fc2+fc'2)r2/8 

where < h^h^hh >= |/3|^, 

Q = q- rdt'F[r(q,t')],Q' = q'- T rft'F[r(q', t'), 

JO JO 

and the quantity r2 is defined (similarly to ri) as r| = r^/ll + Ar'^X'^'^). 

The rest of the calculations required for obtaining <: P :> can be carried 
out according to the scheme outlined for the case of < / >. 

Up to now, we have dealt with the equal-time correlation function <: I{t)I{t) :>. 
This analysis is relevant to the experimental situation in which the response 
time of the detector, t^, is much less than the characteristic time of the phase 
diffuser, r^. In what follows, we consider the opposite case, r^ >> r^. The 
detector averages intensity fluctuations during the time interval r^. Lower- 
ing of the noise level may be expected for this detector. Although the time 
interval, r^, is much larger than, r^, at the same time it should be much 
shorter than the characteristic time of the turbulence evolution (frozen tur- 
bulence), i.e. Td « Ta = l/va, where I and Va are the characteristic radius 
of the turbulent eddies and their transverse flow velocity across the beam, 
respectively. The average value <: I{t)I(t + r) :> with t r^ r^ gives the 
dominant contribution to the quantity measured by this detector. In this 
case, the atmospheric conditions of light propagations may be considered to 
be fixed, while the initial correlations of four field operators 6q,6q, which 
enter <: I{t)I{t + r) :>, should be calculated accounting for the very dif- 
ferent random phases (^(to) and ipito + r) introduced by the diffuser. These 
phases do not correlate (< ip{Y,t)ip{Y\t + r) >= 0) even at r = r'. Hence, 
averaging over the random phases of each of the product of four operators is 
reduced to calculations of < e-^{^['-i(*o)]+^['-2(to+r)]-v'[r3(to+r)]-^[r4(to)]} > ^haX 
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is equal to < e-^^^('"i)-'^('"'')> >< e-^^'^(''2)-'^(''«)> >. As a result, the general 
expression for the intensity correlations will become different from Eq. (|3Up. 
The difference is that Eq. (1311) should become 




bb > { e"(^ +^ )^i/2-(^ +^ )''o/8 (32) 



_^g-[(Q-Q')' + (k+k')V4]rg/4-[(Q+Q')2 + (k-k')V4]r2/4 



The two terms in the braces of Eq. ()32|) correspond to the two summands 
in Eq. ()29|) which describe two types of wave correlations in the course 
of four waves propagation through turbulent atmosphere. As we see, the 
contributions of both trajectories are not equivalent when ri 7^ tq. The 
effect is entirely due to partial coherence and may be controlled by means of 
variation of ri. For the case of ri = r^, Eq. (jH^ is reduced to Eq. (jSH). 

The relative contributions of the two terms in Eq. P^ to the intensity 
correlation function can be easily estimated. The effective volumes of inte- 
gration of the first and second terms over q and q' are of the order of rf 
and t^'^Tq'^, respectively. Moreover, for long-distance propagation, the inte- 
grations over k and k' are confined to the "turbulent" terms, but not by rf ^ 
or r^^. Hence the second term gives a contribution, which is (ro/ri)^ times 
less than the first one. It follows from a comparison of Eqs. ()31|) and (jH^ 
that the intensity correlation function measured by a fast detector is approx- 
imately twice the value of a slow-detector measurement when tq >> ri. Of 
course, this estimate is only valid for long-distance propagation. The next 
section deals with the case of a slow detector in more detail. 



4 Calculations of the intensity fluctuations 

It follows from previous considerations that the intensity correlation function 
for a slow detector is given by 

<: lit + T)I{t) :>= f ^""fyn < 6+6+66 > x (33) 
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^g-i{k[r-c(q)t]+k'[r-c(q')t] + (c/9o) Jo di'i' (kF[r(q,J')]+k'F[r(q',J')])} > ^ 

As we see, the problem is reduced to averaging over the refractive index 
fluctuations and many-fold integrations. It is worthwhile to recall here that 
Q and Q' are dependent on the fluctuating force F. To get a linear form 
with respect to F in the exponents of Eq. (jHHj) . the integral representations 
for the factors containing Q and Q' may be used: 



g-Q2r2/2 ^ / ^^gipQ-pV(2r-2)^ (■g^^ 



1 



^_(Q+Q')2r2/4-(Q-Q')2r2/4 ^ f ^P^P^gipQ+ip'Q'-(p-p')V(4rg)-{p+p')V(4r2)_ 

J (27rrori)2 

(35) 
After substitution of expressions (jH^ and (jHSI) into Eq. (jSHl), further anal- 
ysis is facilitated considerably. In this case the fluctuating field enters the 
intensity correlation function only via the factor 

^_iJ-Jdt'{(p+kt's/go)F[r(q,t')] + (p'+k't's/<?o)F[r(q',t')]} r^Q\ 

The averaging of Eq. (j36p over the refractive index fluctuations may be per- 
formed straightforwardly if we consider the trajectories r(q, t') and r(q',t') 
to be unperturbed by the turbulence, i.e. r(q, t') = r + c(q)(t' — t), and 
r(q', t') = r + c(q')(t' — t). Then, the average value of Eq. (IHUj) is given by 



0/ 



expl - OmSCyql r dx H dgg-^'^e-^'"''/^^^' [(p + \ix/qof + (p + k'x/g, 

(37) 
+2(pi|+fci|a:/go)(pj| + fc||a;/go)(^o-^2)+2(p±+/cxx/go)(Px + ^Wgo)(^o+^2)]|, 

where Jq and J2 are zeroth and second order Bessel functions with the ar- 
guments equal to g\(\ — q'|(-2 — x)/qo. The indices {y} and {_l} indicate the 
parallel and perpendicular to q — q' components of the corresponding 2D 
vectors. In the derivation of Eq. (j37|) we have used the relations 

< F[r(q,t')]F[r(q,0] >=< F[r(q, t')-r(q, t")]F(0) >=< F[c^{t'-t")]F{0) > 

(38) 
and 

< F[r(q, t')]F[r(q', t")] >=< F[r(q, t') - r(q', t")]F(0) > (39) 
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=< F[cq,(t'-t")-Cq-q'(i-i')]F(0) >^< F[ce,{t' -f") - c^^^^t -t')]F{0) >, 

where e^ is a unit vector in the z-direction. The two first terms in square 
brackets of Eq. ()37p describe correlations of waves with the same q and 
q', while two other terms describe cross-correlations. The contribution of 
the last terms decreases with increasing distance of propagation, z, because 
<^,2(z/) -^ when y ^ oo. Neglecting the cross-correlation terms, we can 
easily obtain the asymptotic value {z -^ oo) of the intensity fluctuations and 
the scintillation index. The scintillation index is given by 

< / >2 rf{ 



Similar result for the limiting case — ^ was obtained in ^8]. At the 
same time, this result differs from the asymptotic value equal to 1, obtained 
in [21]. The difference may arise because of a different assumptions on the 
coherent properties of the sourse used by the authors of [21] . Namely, in [21] 
the authors consider a small coherence time equal to inverse bandwidth of 
the laser generation, while our results are valid for small characteristic times 
of local phase fluctuations introduced by the dynamic phase screen. 

As we see, the scintillation index tends to zero when z — > oo, (ri/ro)^ — *• 
0. This property of partially coherent radiation is favorable for practical 
utilization. 

Similar reasonings may be used to obtain the correlator of arbitrary (the 
r/ith) order. It is given by 







(41) 



where Tg < \ti — tj\ < t^ and am = 1,2,9,44,265.. for m = 2,3,4,5,6.., 
respectively. 

Taking into account Eqs. ([Hn jl - flHU]) and all terms in Eq. (|H7|) we obtain an 
analytic expression for the intensity fluctuations in which many integrations 
can be performed analytically. The rest (three-fold integral) can be evaluated 
numerically. This allows us to analyze the effect of turbulence on a^ at large 
but flnite distances. 

Eq. (jH7|) is derived under the assumption that the trajectories r(q, t') 
in the fluctuating force F are straight lines. Let us analyse the effect of 
distortion of trajectories due to the fluctuating force. With regard for the fact 
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that the force F is responsible for changing the photon transverse momentum 
(and the transverse velocity), the account for the dependence of F on F means 
the account of nonlinear in F effects in the photon transverse displacements. 
These effects are beyond the applicability of the fourth-moment equation pP , 
0, therefore, their studies are impossible by means of the Yakushkin method. 
Also, this means that papers ^H], QH]; which are based on the Yakushkin 
method, stay this problem out of consideration. 

First of all, it should be noted that the correlation function < F(ri)F(r2) > 
depends on ri — r2. If both points belong to the same trajectory [ri = 
r(q, t'),r2 = r(q, t")] this difference is equal to r(q, t') — r(q, t") ^ cez{t' — t"). 
[See Eq. (IHH|).] Consequently \t' — t"\ < Lq/c and only negligible particle dis- 
placements in perpendicular to the z-direction may occur for such very short 
time intervals. In the other case, when two points belong to different trajec- 
tories [ri = r(q, t'),r2 = r(q',t")] we have 



ri-r2 ^ce,(t'-t")-Cq-q'(t-t') + - / rfti(ti-t'){F[r(q, ti)]-F[r(q', ti)]}. 

qoJt' 

(42) 
It can be easily seen from Eq. ()42|) that for any given \t' — t"\ < Lq/c, 
the last two terms on the right-hand side may be comparable with the first 
term for sufficiently large values of t — t'. In this case, the deviations of the 
trajectories from straight lines, which enter arguments in the left hand side 
of Eq. (jnni), should be taken into account. One can say that the deviations 
are accumulated throughout the whole propagation path and, in contrast to 
the case of the same trajectory, may become sufficiently large to influence 
the wave correlations. 

For further analysis, an important property of the correlation function 
< F(ri — r2)F >, where ri — r2 is given by Eq. (02), has to be noted. There is 
almost no correlation between F entering the argument of F and the function 
itself. This is because of the negligibly small time intervals (or characteristic 
distances) where these functions correlate: < ti — t', t" < Lq/c. Then, if 
we consider as previously that Lq « z, the averaging of < F(ri — r2)F > 
may be undertaken in two steps: flrstly, we average this quantity considering 
ri — r2 as a flxed parameter and after that the remaining averaging should 
be performed. The result is 

< F,(ri - r2)F,(0) >= J dgg,g,tp{g) < e-'^^'^-^'^ >, (43) 

where g is a three-dimensional vector and the indices i,j denote the compo- 
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nents perpendicular to z. The procedure described above may be repeated 
to obtain < e"*^^'"^"'"^^ >. Again, this function may be expressed in terms 
of two-point correlation functions in which the points belong to the same 
or different trajectories. At this stage we will simplify further analysis by 
imposing the approximation 

< e-^^("^^''2) >^< e-*s"^ >< e'^^' >, (44) 

where the correlation of different trajectories is neglected. Physically, this 
approximation requires that two photons with different transverse momenta 
propagate in spatial areas with different refractive indices and experience 
different fluctuating forces. Therefore they gain different values of the trans- 
verse velocity and transverse displacements. It is evident that for any finite 
value of q — q', the displacements become uncorrelated when z ^* oo. 

Within the approximation (j44j) . the effect of refractive index fluctuations 
on trajectories in Eq. (jHUj) may be accounted for by multiplying all the Bessel 
functions in Eq. (jHTj) by the factor e~"'^^~^' , where 

7 = 7(^) = 0.011vr^(2vr)V3r(^i)c^/-V3^^ 

This factor reduces the effect of correlations of different trajectories by taking 
into account randomization of the particle displacements from straight lines. 



5 Weak irradiance fluctuations 

For the case of weak turbulence or short distance propagation, the beam 
characteristics are the same as for propagation in free space. Small devia- 
tions from free space regime may be accounted using perturbation methods. 
As previously, the intensity fluctuations may be described in terms of the 
function (^(q, q'; k, k'). Further analysis is simplified if one uses an iterative 
procedure not for G, but for its fluctuating part, T, defined as 

T(q, q'; k, k') = ^(q, q'; kk')- < 6q+k/2Vk/2 >< &J+k'/2V-k'/2 > • (45) 
The equation of motion for T is given by 
[9,-zA(q,q';k,k')]T(q,q';k,k')=^a;oE^g{'^(q + g/2,q';k-g,k') (46) 
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-T(q - g/2, q'; k - g, k') + [q, k ^ q', k']}, 

where the symbol [q, k <->^ q', k'] indicates two summands, which are similar 
to previous ones but with the variables interchanged as indicated, and 

A(q, q'; k, k') = C^q+k/2 + ^q'+k'/2 — ^q'-k'/2 — ^q-k/2- 

It follows from Eq. pUj) that the average value of T may be written as 

< T(q, q'; k, k') >= luJo V /* rft'e*A(q,q';k,k'){t-t') ^47) 

X < ng{T(q+g/2, q'; k-g, k')-T(q-g/2, q'; k-g, k') + [q, k ^ q', k']}|,=,, >, 

where the initial condition < T(q, q';k, k') > |j=o = corresponding to the 
case of slow detector, was used. The right-hand side of Eq. (jTTj) . dependent 
on average value of < n^T >, is still unknown. An expression similar to Eq. 
(I47|) may be derived again for < n^T >. Then, we will simplify the problem 
by considering the time dependence of all operators in the integrand to be 
determined by free-space propagation laws. The approach is in the spirit of 
the classical Rytov approximation IJ. This will make it possible to integrate 
all terms and to get explicit forms for < rZgT > as well as for < T >. Using 
some simple algebra, it can be proved that the scintillation index is given by 

a' = alK{z,ro,r,), (48) 

where af is the Rytov variance defined by af = 1.23C^qQ z^^^^ and 



K{z,po,pi) = 4.24 dr dxx ^'^expl 



x^ 



%^ , ^2 Po + pI 

+ T 



,2^2 



A-n^z 4 + pgpf 



(49) 



xsin —— 



V 2 A+pipir 

where pg ^ = r^^qQ/z. As done previously, in course of derivation of Eq. 
(|48p . we have considered the propagation distance to be much longer than 
the characteristic scale of turbulence variation [z ^ Lq). It follows from 
Eq. (jlUj) that in the limit of pi -^ 00, Iq -^ 0, we have the result of Rytov 
theory (a^ = af) because K -^ 1. The quantity in the square brackets of 
Eq. ()49|) is negligible in this case. With decreasing initial beam coherence, 
Vi becomes smaller and it may occur that the quantity in square brackets 
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Figure 2: Dependence of scintillation index on turbulence strength for weak 
fluctuations; go and Iq are as in Fig. 1. Decreasing of scintillation index with 
the decrease of source coherence is seen (similar to results of O. Korotkova 
et al. ini)- 



becomes sufficiently large to influence the result of integration. The character 
of wave propagation becomes modified from the plane wave regime to the 
spherical wave regime. This is accompanied by a decrease of o"^. The effect 
is saturated at some small values of ri ~ Izlr^q^ and further decrease of the 
coherent length Ac has no effect on a^. Also, when po is small, the effect of r\ 
variation is of no significance. Therefore, there is the opportunity to control 
cr^ but this is only possible at a sufficiently large aperture radius tq and small 
distance z. Furthermore, the reduction of o"^ is limited by some finite value. 
Fig. 2 illustrates the effect of partial coherence just for the favorable case 
when Po = 5. 

6 Discussion 

Figs. 3 and 4 show the dependence of scintillation index on the turbulence for 
two different distances, z. A well-pronounced effect for the decreasing of o"^ 
with the decrease of the initial coherence can be seen in the range of strong 
turbulence. There is a very simple physical explanation of the reduction 
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Figure 3: Plots o"^(C^) for different initial coherence (ri/ro)^; go and /q are 
as in Fig. 1. Left-side curves calculated with employing Eqs. ()48p and 
P^ exhibit negligible effect of partial coherence (merged curves), that is in 
contrast to results shown in Fig. 2. Two very different parameters pg are 
used in these cases. Black squares show o^ obtained with the assumption of 
straight trajectories in Eq. (|^. 



of cr^. Two things are very important for understanding this phenomenon. 
First of all, in the course of the irradiance propagation the beam acquires 
the properties of Gaussian statistics. Therefore, the asymptotic value of the 
intensity correlations < I{t)I{t) > is given by 2 < / >^. On the other 
hand, the quadratic detector counts are determined not by simultaneous 
correlations, but by the average < I{t)I{t + r) >, where the characteristic 
value of r is of the order of r^. In the limiting case of Ac ^ and r » Tg 
(slow detector), there are no correlations between I{t) and I{t + r). Hence, 
< I{t)I{t + r) >=< I >^, and the normalized variance of the intensity 
fluctuations is negligible (cr^ -^ 0). This physical picture is quite similar 
to the well-known Hanbary-Braun-Twiss effect [22] 5 or photon bunching for 
thermal light: at zero delay the correlation function has twice the value for 
long delays. For finite values of Ac, the irradiance differs from thermal light 
resulting in a finite value of the scintillation index. Our theory gives the 
scintillation index equal to r^/rg, that is the relative part of the aperture 
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Figure 4: The same as in Fig. 3 but for longer distance. The upper curve 
corresponds to coherent beam. It approaches asymptotically to value cr^ = 1. 



where the exiting light may be considered as a coherent light. 

It should be noted that the finite value of a^ shows an absence of full 
thermalization because the memory of the photon fiux about source correla- 
tions still exists. This memory may be lost due to fiuctuations of the photon 
transit time 6t. [The transit time was assumed to be constant (equal to z/c) 
in the previous analysis.] When 6t ^ 0, the initial conditions for /(r, q, t) 
correspond not to its value at the aperture plane, but, for example, to some 
outlying point Sz = c\6t\ where the turbulence has already modified the 
waves. The estimate of 6z is given by its rms value. Setting 



Sz 



we can easily obtain < (Sz)'^ >^ 
Lo 



/ dz5n{z) 
Jo 



0M67t'^zCIlI^\ 



When C^ = 10' 



13 



m 



-2/3 



lO^m, we have < {5z)'^ >V2^ o.55 * 10"^m and \5t\ 



2 * 10 ^^s, which is negligibly small. 



10m, z 

Other causes of transit time fiuctuations arise from fiuctuations of the 
transverse velocity of the photons. For nonzero q± the photon velocity in the 
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z direction is given by 

Assuming g_|_ to be caused by the turbulence, we have 

q(z)= rF{z')dz'/c, (51) 

Jo 

where indices _l are omitted again. Then, using Eqs. (jSUjl and ()51|) we may 
estimate 6t as 

For the same parameters used previously, (5t ~ 3 * lO^^^s and < {Sz^ >^/^~ 
9 * 10~^?7i, which is negligible again. 

The most serious assumption of our approach is that the influence of the 
turbulence on the photon distribution can be described in terms of a random 
force F(r) (see Eq. ((Tj)) modifying photon momentum q±2;. Mathematically, 
this approximation can be justified when the width of the photon distribution 
in momentum space is greater than those turbulence wave vectors, which give 
the dominant contribution to cr^. In the vicinity of the source, the photon 
momentum is distributed within the range of the order of vr/ri. (See Eq. 
(J2S1)-) The turbulence spectrum covers very broad interval of wave vectors 
[tt/Lq, 7r//o] where Lq/Iq ~ 10^. Therefore, it is not clear a priori what wave 
vector should be taken for the comparison. Of course, the condition ri < Iq is 
sufficient to validate our approach. At the same time, this condition imposes 
very rigid restriction on real optical systems. Fortunately, it is not obligatory 
for a reliable solution for long-distance propagation. Fist of all, it may be 
assumed that a broad range of the turbulence spectrum, rather than wave 
vectors equal to vt/Zq, contributes significantly to measured quantities. If this 
is true, these quantities would not be sensitive to the boundary values of wave 
vectors. In this context, it is important to note that the beam radius really 
does exhibit a weak dependence on Iq. [It depends on /q through T r^ Iq , 
see Eq. fl26|) .] Hence, it is reasonable to consider the characteristic value of 
the turbulence wave vector to be much smaller than tt/Zq- 

On the other hand, due to the action of the random force, there is a 
diffusion-like increase in the transverse momentum of photons in the course 
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of their propagation. It is just that quantity (not 7r/ri) which should be 
compared with the turbulence wave vectors for large z. Using Eq. (|^ we 
may estimate the increase of the transverse momentum as 

< g2 >= 0.0667r2r('i')go'^C^(27r//o)^/^ (53) 

Substituting the previous parameters into Eq. (j3^ . we get < (f > /tt^/q"^ ~ 
10^. This estimate shows our approach to be reliable for long distance prop- 
agation. At small distances, the perturbation theory is applicable. 

7 Conclusion 

The approach presented in this paper can be used for both stationary beams 
as well as for beams with varying intensity. Also, the statistics of the exit- 
ing irradiance may differ from the statistics of coherent or partially coherent 
beams. Our method of obtaining the correlation function for the intensity 
fluctuations is not based on the solution of an equation for the fourth-order 
correlation function. Moreover, we have obtained the asymptotic value for all 
moments of the intensity fluctuations. Our method does not employ Markov 
approximation used previously to derive the equation for the fourth-order 
correlation function. Also, we do not use the so-called quadratic approxima- 
tion, which like the Markov approximation is nothing more than an artificial 
modification of the refractive index covariance in order to simplify the theo- 
retical treatment. Both simplifications are inherent to papers ^^ and P^ . 
Besides that, there was assumed /q = 0, that is impossible within our consid- 
eration. Suffice it to say that in the case of /q = 0, the quantities T and 7 lose 
their physical sense because of infinitely large values. In spite of the presence 
of so serious approximations, the asymptotic (2; — ;> cx) or C^ — ;> 00) results of 
[TH] and our paper results are very similar. The explanation of this follows 
from the result obtained: the asymptotic value of a^ does not depend on the 
atmosphere turbulence mechanism. Therefore, it does not matter what kind 
of turbulence spectrum is employed in the case of saturated fluctuations. At 
the same time, our calculations show (see Figs. 3 and 4) that the distinct 
asymptotic regime scarcely can be reached in real experiments. Thus, for 
practice, it is more important to know the behavior of a^ at large but flnite 
C^ and z than its asymptotic value. In this case, the real turbulence spec- 
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trum as well as an adequate theoretical analysis are required. In this context, 
the motivations of our studies are evident. 

Our calculations show how to suppress the intensity fluctuations using 
a PCB. The conditions required for this are a random phase modulation 
of the signal and the use of a slow detector. Estimations of the utility of 
the PCB should take into account the negative factor of worsening for the 
receiving system resolution when t^ is increased. Also, it is important for 
the experiment performance, that the ratio Tg/rd is always finite, and just 
this ratio rather than (ri/ro)^ determines the asymptotic value of a when 
ijs/T~d) > {ri/roY. Hence, the regime of fast phase modulations (but not the 
increase of t^) is preferable in both cases. 
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